## Warning: st_centroid assumes attributes are constant over geometries
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in validateCoords(lng, lat, funcName): Data contains 1 rows with either
## missing or invalid lat/lon values and will be ignored
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
Below are mapview interactive maps of raw data
mapview(grid_110m , alpha.regions = 0.2, alpha = 1 , col.regions="darkgrey" ) + mapview(grid_22m , alpha.regions = 0.2, alpha = 1 , col.regions="grey") + mapview(focals_data_cropped , cex=2, zcol="nutcrackin") + mapview(nutcracking_data_cropped, col.regions ="violet", cex=2 ) + mapview(palm_trees, col.regions = "green3", cex=2)
mapview(focals_data_cropped , cex=2, zcol="individual")
grid_110m
## Simple feature collection with 1000 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 288756.2 ymin: 9576576 xmax: 290076.2 ymax: 9577126
## Projected CRS: WGS 84 / UTM zone 24S
## First 10 features:
## grid_id id left top right bottom cell_id count_palm_trees
## 1 1 1 289966.2 9576906 290076.2 9576796 483 0
## 2 1 1 289966.2 9576906 290076.2 9576796 601 0
## 3 1 1 289966.2 9576906 290076.2 9576796 482 2
## 4 1 1 289966.2 9576906 290076.2 9576796 484 0
## 5 1 1 289966.2 9576906 290076.2 9576796 602 0
## 6 1 1 289966.2 9576906 290076.2 9576796 603 0
## 7 1 1 289966.2 9576906 290076.2 9576796 604 0
## 8 1 1 289966.2 9576906 290076.2 9576796 542 1
## 9 1 1 289966.2 9576906 290076.2 9576796 543 0
## 10 1 1 289966.2 9576906 290076.2 9576796 544 0
## rawstones_average geometry
## 1 1.33 POLYGON ((289966.2 9576906,...
## 2 1.33 POLYGON ((289966.2 9576906,...
## 3 1.33 POLYGON ((289966.2 9576906,...
## 4 1.33 POLYGON ((289966.2 9576906,...
## 5 1.33 POLYGON ((289966.2 9576906,...
## 6 1.33 POLYGON ((289966.2 9576906,...
## 7 1.33 POLYGON ((289966.2 9576906,...
## 8 1.33 POLYGON ((289966.2 9576906,...
## 9 1.33 POLYGON ((289966.2 9576906,...
## 10 1.33 POLYGON ((289966.2 9576906,...
mynewspdf <- merge(grid_110m,average_stones_plot, duplicateGeoms = T)
For this, I will one model out of many I ran, that is the most complete and drives home the gist. This gives similar results to one that accounts for the length of follow (i.e. 1 to 5 minutes), but we can wait till the the actual tool use rates per follow are calculated. Essentially, this is a logistic GLMM that models the probability of using tools in a follow as predicted by empirical palm tree density in the 22m cell where the follow started (double check), the posterior density of mean stone density in the corresponding 110m cell, and an interaction between the 2.
We estimated varying intercepts and slopes for all parameters for each individual capuchin. Stone density was models as a Poisson GLMM with a log link, where we estimated a unique varying intercept for each 110m grid cell. This estimate (or a centered version of it) gets passed to the main GLMM, so we can use the distribution of stones to inform the mean, rather than a single, overconfident point estimate.
set.seed(345)
m10me <- ulam(
alist(
#stone model
stone_raw_grid ~ dpois( lambda ), #define poisson likelihood
log(lambda) <- a_bar + a_grid[stone_grid_index],
#priors
a_bar ~ dnorm( 3 , 2 ), #mean effect on logscale
a_grid[stone_grid_index] ~ dnorm( 0 , sigma_stone ), #priovarying effects centered on mean
sigma_stone ~ dexp(1), #p
#nutcracking model
nutcrackin ~ dbinom(1, p ),
logit(p) <-ap + a_id[id,1] + (bS + a_id[id,2])*a_grid[grid_id_follow] +
(bPT + a_id[id,3])*count_palm_std + (bPTxS + a_id[id,4])*a_grid[grid_id_follow]*count_palm_std,
# adaptive priors - non-centered
transpars> matrix[id,4]:a_id <-
compose_noncentered( sigma_id , L_Rho_id , z_id ),
matrix[4,id]:z_id ~ normal( 0 , 1 ),
# fixed priors
c(ap,bS,bPT,bPTxS) ~ dnorm( 0 , 1 ),
vector[4]:sigma_id ~ dexp(1),
cholesky_factor_corr[4]:L_Rho_id ~ lkj_corr_cholesky( 3 ),
# compute ordinary correlation matrixes from Cholesky factors
gq> matrix[4,4]:Rho_id <<- Chol_to_Corr(L_Rho_id)
) , data=data_list_daily , chains=4 , cores=4, log_lik=FALSE, iter=4000, pars_omit = 'z_id' , control=list(adapt_delta=0.99))
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 2 Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 3 Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 4 Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[3] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f547f39a7f8.stan', line 36, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 2 Iteration: 100 / 4000 [ 2%] (Warmup)
## Chain 3 Iteration: 100 / 4000 [ 2%] (Warmup)
## Chain 1 Iteration: 100 / 4000 [ 2%] (Warmup)
## Chain 4 Iteration: 100 / 4000 [ 2%] (Warmup)
## Chain 3 Iteration: 200 / 4000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 4000 [ 5%] (Warmup)
## Chain 3 Iteration: 300 / 4000 [ 7%] (Warmup)
## Chain 4 Iteration: 200 / 4000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 4000 [ 5%] (Warmup)
## Chain 2 Iteration: 300 / 4000 [ 7%] (Warmup)
## Chain 3 Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 4000 [ 7%] (Warmup)
## Chain 3 Iteration: 500 / 4000 [ 12%] (Warmup)
## Chain 2 Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4 Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3 Iteration: 600 / 4000 [ 15%] (Warmup)
## Chain 1 Iteration: 300 / 4000 [ 7%] (Warmup)
## Chain 3 Iteration: 700 / 4000 [ 17%] (Warmup)
## Chain 2 Iteration: 500 / 4000 [ 12%] (Warmup)
## Chain 4 Iteration: 500 / 4000 [ 12%] (Warmup)
## Chain 1 Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4 Iteration: 600 / 4000 [ 15%] (Warmup)
## Chain 2 Iteration: 600 / 4000 [ 15%] (Warmup)
## Chain 3 Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4 Iteration: 700 / 4000 [ 17%] (Warmup)
## Chain 2 Iteration: 700 / 4000 [ 17%] (Warmup)
## Chain 4 Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3 Iteration: 900 / 4000 [ 22%] (Warmup)
## Chain 2 Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4 Iteration: 900 / 4000 [ 22%] (Warmup)
## Chain 1 Iteration: 500 / 4000 [ 12%] (Warmup)
## Chain 4 Iteration: 1000 / 4000 [ 25%] (Warmup)
## Chain 3 Iteration: 1000 / 4000 [ 25%] (Warmup)
## Chain 2 Iteration: 900 / 4000 [ 22%] (Warmup)
## Chain 4 Iteration: 1100 / 4000 [ 27%] (Warmup)
## Chain 1 Iteration: 600 / 4000 [ 15%] (Warmup)
## Chain 2 Iteration: 1000 / 4000 [ 25%] (Warmup)
## Chain 4 Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3 Iteration: 1100 / 4000 [ 27%] (Warmup)
## Chain 2 Iteration: 1100 / 4000 [ 27%] (Warmup)
## Chain 1 Iteration: 700 / 4000 [ 17%] (Warmup)
## Chain 4 Iteration: 1300 / 4000 [ 32%] (Warmup)
## Chain 2 Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4 Iteration: 1400 / 4000 [ 35%] (Warmup)
## Chain 3 Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1 Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2 Iteration: 1300 / 4000 [ 32%] (Warmup)
## Chain 4 Iteration: 1500 / 4000 [ 37%] (Warmup)
## Chain 2 Iteration: 1400 / 4000 [ 35%] (Warmup)
## Chain 3 Iteration: 1300 / 4000 [ 32%] (Warmup)
## Chain 4 Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4 Iteration: 1700 / 4000 [ 42%] (Warmup)
## Chain 2 Iteration: 1500 / 4000 [ 37%] (Warmup)
## Chain 1 Iteration: 900 / 4000 [ 22%] (Warmup)
## Chain 3 Iteration: 1400 / 4000 [ 35%] (Warmup)
## Chain 4 Iteration: 1800 / 4000 [ 45%] (Warmup)
## Chain 2 Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4 Iteration: 1900 / 4000 [ 47%] (Warmup)
## Chain 1 Iteration: 1000 / 4000 [ 25%] (Warmup)
## Chain 3 Iteration: 1500 / 4000 [ 37%] (Warmup)
## Chain 2 Iteration: 1700 / 4000 [ 42%] (Warmup)
## Chain 4 Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4 Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 4000 [ 27%] (Warmup)
## Chain 3 Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2 Iteration: 1800 / 4000 [ 45%] (Warmup)
## Chain 4 Iteration: 2100 / 4000 [ 52%] (Sampling)
## Chain 1 Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2 Iteration: 1900 / 4000 [ 47%] (Warmup)
## Chain 3 Iteration: 1700 / 4000 [ 42%] (Warmup)
## Chain 1 Iteration: 1300 / 4000 [ 32%] (Warmup)
## Chain 4 Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2 Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2 Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3 Iteration: 1800 / 4000 [ 45%] (Warmup)
## Chain 1 Iteration: 1400 / 4000 [ 35%] (Warmup)
## Chain 4 Iteration: 2300 / 4000 [ 57%] (Sampling)
## Chain 2 Iteration: 2100 / 4000 [ 52%] (Sampling)
## Chain 3 Iteration: 1900 / 4000 [ 47%] (Warmup)
## Chain 1 Iteration: 1500 / 4000 [ 37%] (Warmup)
## Chain 4 Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2 Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1 Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4 Iteration: 2500 / 4000 [ 62%] (Sampling)
## Chain 2 Iteration: 2300 / 4000 [ 57%] (Sampling)
## Chain 3 Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3 Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1 Iteration: 1700 / 4000 [ 42%] (Warmup)
## Chain 4 Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2 Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1 Iteration: 1800 / 4000 [ 45%] (Warmup)
## Chain 3 Iteration: 2100 / 4000 [ 52%] (Sampling)
## Chain 4 Iteration: 2700 / 4000 [ 67%] (Sampling)
## Chain 2 Iteration: 2500 / 4000 [ 62%] (Sampling)
## Chain 1 Iteration: 1900 / 4000 [ 47%] (Warmup)
## Chain 4 Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3 Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2 Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 4 Iteration: 2900 / 4000 [ 72%] (Sampling)
## Chain 1 Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1 Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2 Iteration: 2700 / 4000 [ 67%] (Sampling)
## Chain 3 Iteration: 2300 / 4000 [ 57%] (Sampling)
## Chain 4 Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1 Iteration: 2100 / 4000 [ 52%] (Sampling)
## Chain 2 Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4 Iteration: 3100 / 4000 [ 77%] (Sampling)
## Chain 3 Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1 Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2 Iteration: 2900 / 4000 [ 72%] (Sampling)
## Chain 4 Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3 Iteration: 2500 / 4000 [ 62%] (Sampling)
## Chain 1 Iteration: 2300 / 4000 [ 57%] (Sampling)
## Chain 2 Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 4 Iteration: 3300 / 4000 [ 82%] (Sampling)
## Chain 2 Iteration: 3100 / 4000 [ 77%] (Sampling)
## Chain 3 Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1 Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4 Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2 Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3 Iteration: 2700 / 4000 [ 67%] (Sampling)
## Chain 1 Iteration: 2500 / 4000 [ 62%] (Sampling)
## Chain 4 Iteration: 3500 / 4000 [ 87%] (Sampling)
## Chain 2 Iteration: 3300 / 4000 [ 82%] (Sampling)
## Chain 3 Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1 Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 4 Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2 Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1 Iteration: 2700 / 4000 [ 67%] (Sampling)
## Chain 3 Iteration: 2900 / 4000 [ 72%] (Sampling)
## Chain 4 Iteration: 3700 / 4000 [ 92%] (Sampling)
## Chain 2 Iteration: 3500 / 4000 [ 87%] (Sampling)
## Chain 1 Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4 Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 3 Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2 Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1 Iteration: 2900 / 4000 [ 72%] (Sampling)
## Chain 4 Iteration: 3900 / 4000 [ 97%] (Sampling)
## Chain 2 Iteration: 3700 / 4000 [ 92%] (Sampling)
## Chain 3 Iteration: 3100 / 4000 [ 77%] (Sampling)
## Chain 1 Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 4 Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4 finished in 237.3 seconds.
## Chain 2 Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 3 Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1 Iteration: 3100 / 4000 [ 77%] (Sampling)
## Chain 2 Iteration: 3900 / 4000 [ 97%] (Sampling)
## Chain 1 Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3 Iteration: 3300 / 4000 [ 82%] (Sampling)
## Chain 2 Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2 finished in 252.6 seconds.
## Chain 1 Iteration: 3300 / 4000 [ 82%] (Sampling)
## Chain 3 Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1 Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 3 Iteration: 3500 / 4000 [ 87%] (Sampling)
## Chain 1 Iteration: 3500 / 4000 [ 87%] (Sampling)
## Chain 3 Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1 Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3 Iteration: 3700 / 4000 [ 92%] (Sampling)
## Chain 1 Iteration: 3700 / 4000 [ 92%] (Sampling)
## Chain 3 Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1 Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 3 Iteration: 3900 / 4000 [ 97%] (Sampling)
## Chain 1 Iteration: 3900 / 4000 [ 97%] (Sampling)
## Chain 3 Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3 finished in 300.5 seconds.
## Chain 1 Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1 finished in 300.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 272.8 seconds.
## Total execution time: 301.1 seconds.
Lets look at the parameter outputs. First the main effects.
#precis(m10me)
plot(precis(m10me))
## 144 vector or matrix parameters hidden. Use depth=2 to show them.
#hetergenity amongs grids in stones
plot(precis(m10me , depth=3 , pars="a_grid"))
Tree density has a reliable positive impact on the probability of using
tools, whereas stone abundance has a small, unreliable positive impact
on tool use. There is a weak interaction effect as well- places where
there are more .
post <- extract.samples(m10me)
str(post)
## List of 12
## $ a_bar : num [1:8000, 1] 2.47 2.47 2.52 2.53 2.68 ...
## $ a_grid : num [1:8000, 1:40] -1.69 -1.79 -1.39 -1.65 -2.55 ...
## $ sigma_stone: num [1:8000, 1] 1.33 1.14 1.15 1.23 1.43 ...
## $ z_id : num [1:8000, 1:4, 1:17] 0.154 -1.301 -1.353 -1.242 -0.328 ...
## $ bPTxS : num [1:8000, 1] -0.00251 -0.00142 -0.06499 0.02509 0.07459 ...
## $ bPT : num [1:8000, 1] 0.606 0.773 0.748 0.696 0.661 ...
## $ bS : num [1:8000, 1] 0.018 0.096 0.1518 0.1415 -0.0304 ...
## $ ap : num [1:8000, 1] -1.33 -1.24 -1.24 -1.17 -1.37 ...
## $ sigma_id : num [1:8000, 1:4] 0.816 0.47 0.477 0.524 0.735 ...
## $ L_Rho_id : num [1:8000, 1:4, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
## $ a_id : num [1:8000, 1:17, 1:4] 0.126 -0.612 -0.645 -0.651 -0.241 ...
## $ Rho_id : num [1:8000, 1:4, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "source")= chr "ulam posterior from object"
p_link_palm <- function(x){
scapula <- mean(post$a_grid)
id <- 1
logit_scale <- with(post , ap + a_id[,id,1]*0 + (bS + a_id[,id,2])*scapula +
(bPT + a_id[,id,3]*0)*x + (bPTxS + a_id[,id,4]*0)*scapula*x )
return(logistic(logit_scale))
}
palm_seq <- seq(from=min(data_list_daily$count_palm_std) ,
to=max(data_list_daily$count_palm_std) , length=30)
p_raw <- sapply( palm_seq, function(x) p_link_palm(x) )
plot(data_list_daily$nutcrackin ~ data_list_daily$count_palm_std ,
pch=19 , col=col.alpha("green4",0.004), ylab="probability of nut cracking",
xlab= "palm tree density/ 22 m^2 grid (standardized)")
pred_median <- apply(p_raw , 2 , median)
lines(pred_median ~ palm_seq , lw=2, col=1 , lty=1)
for (j in sample( c(1:2000) , 100) ){
lines( p_raw[j,] ~ palm_seq , lw=3, col=col.alpha("green4", alpha=0.1) , lty=1)
}
We can also plot predictions of varying effects for all 17 individuals who used tools. One of them does not like using tools.
## per id plot
plot(data_list_daily$nutcrackin ~ data_list_daily$count_palm_std ,
pch=19 , col=col.alpha("green4",0.004), ylab="probability of nut cracking",
xlab= "palm tree density/ 22 m^2 grid (standardized)")
#lines(pred_median ~ palm_seq , lw=2, col=1 , lty=1)
for ( i in 1:max(data_list_daily$id)){
dpred <- list(
id=rep(i,30), #list of ones we will replace
#log_avg_stone_std=rep(median(post$a_bar),30),
count_palm_std=palm_seq,
stone_grid_index=rep(1,30),
grid_id_follow=rep(1,30)
)
link2 <- link(m10me, data=dpred , replace=list(stone_grid_index=a_grid_s_z , grid_id_follow=a_grid_f_z) )
pred_median <- apply(link2$p , 2 , median)
lines(pred_median ~ palm_seq , lw=2, col=col.alpha("green4",0.5) , lty=1)
}
post <- extract.samples(m10me)
str(post)
## List of 12
## $ a_bar : num [1:8000, 1] 2.47 2.47 2.52 2.53 2.68 ...
## $ a_grid : num [1:8000, 1:40] -1.69 -1.79 -1.39 -1.65 -2.55 ...
## $ sigma_stone: num [1:8000, 1] 1.33 1.14 1.15 1.23 1.43 ...
## $ z_id : num [1:8000, 1:4, 1:17] 0.154 -1.301 -1.353 -1.242 -0.328 ...
## $ bPTxS : num [1:8000, 1] -0.00251 -0.00142 -0.06499 0.02509 0.07459 ...
## $ bPT : num [1:8000, 1] 0.606 0.773 0.748 0.696 0.661 ...
## $ bS : num [1:8000, 1] 0.018 0.096 0.1518 0.1415 -0.0304 ...
## $ ap : num [1:8000, 1] -1.33 -1.24 -1.24 -1.17 -1.37 ...
## $ sigma_id : num [1:8000, 1:4] 0.816 0.47 0.477 0.524 0.735 ...
## $ L_Rho_id : num [1:8000, 1:4, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
## $ a_id : num [1:8000, 1:17, 1:4] 0.126 -0.612 -0.645 -0.651 -0.241 ...
## $ Rho_id : num [1:8000, 1:4, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "source")= chr "ulam posterior from object"
p_link_stone <- function(x){
id <- 1
logit_scale <- with(post , ap + a_id[,id,1]*0 + (bS + a_id[,id,2])*x +
(bPT + a_id[,id,3]*0)*0 + (bPTxS + a_id[,id,4]*0)*0*x )
return(logistic(logit_scale))
}
stone_seq <- seq(from=min(data_list_daily$log_avg_stone_std) ,
to=max(data_list_daily$log_avg_stone_std) , length=30)
p_raw <- sapply( stone_seq, function(x) p_link_stone(x) )
plot(data_list_daily$nutcrackin ~ data_list_daily$log_avg_stone_std ,
pch=19 , col=col.alpha("darkslateblue",0.004), ylab="probability of nut cracking",
xlab= "log stone density/ 110 m^2 grid (standardized)")
pred_median <- apply(p_raw , 2 , median)
lines(pred_median ~ stone_seq , lw=2, col=1 , lty=1)
for (j in sample( c(1:2000) , 100) ){
lines( p_raw[j,] ~ stone_seq , lw=3, col=col.alpha("darkslateblue", alpha=0.1) , lty=1)
}
There also be a weak interaction, but treat with caution.
##stone abundance
par( mfrow=c(1,3))
stone_samp <- c(-2,0,2)
for(j in 1:3){
plot(data_list_daily$nutcrackin ~ data_list_daily$count_palm_std ,
pch=19 , col=col.alpha("green4",0.004), ylab="probability of nut cracking",
xlab= "palm tree density/ 22 m^2 grid (standardized)" , main=stone_samp[j] )
palm_seq <- seq(from=min(data_list_daily$count_palm_std) ,
to=max(data_list_daily$count_palm_std) , length=30)
dpred <- list(
id=rep(1,30), #list of ones we will replace
log_avg_stone_std=rep(stone_samp[j],30),
count_palm_std=palm_seq,
stone_grid_index=rep(1,30),
grid_id_follow=rep(1,30)
)
link2 <- link(m10me, data=dpred , replace=list(id=a_id_z , stone_grid_index=a_grid_s_z , grid_id_follow=a_grid_f_z) )
pred_median <- apply(link2$p , 2 , median)
lines(pred_median ~ palm_seq , lw=2, col=1 , lty=1)
for (j in sample( c(1:2000) , 100) ){
lines( link2$p[j,] ~ palm_seq , lw=3, col=col.alpha("green4", alpha=0.1) , lty=1)
}
}
#dev.off()
m_h2 <- ulam(
alist(
height_dm ~ dzipois( p , lambda ),
logit(p) <- ap + a_id[id,1] + (bNp + a_id[id,3])*nutcrackin ,
log(lambda) <- al + a_id[id,2] + (bNl + a_id[id,4])*nutcrackin ,
# adaptive priors - non-centered
transpars> matrix[id,4]:a_id <-
compose_noncentered( sigma_id , L_Rho_id , z_id ),
matrix[4,id]:z_id ~ normal( 0 , 1 ),
# fixed priors
c(ap,bNp,bNl) ~ dnorm( 0 , 1 ),
al ~ dnorm( 3.5 , 2 ),
vector[4]:sigma_id ~ dexp(1),
cholesky_factor_corr[4]:L_Rho_id ~ lkj_corr_cholesky( 3 ),
# compute ordinary correlation matrixes from Cholesky factors
gq> matrix[4,4]:Rho_id <<- Chol_to_Corr(L_Rho_id)
) , data=data_list_height , chains=4 , cores=4 , cmdstan = TRUE)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[4] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[4] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[4] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[4] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[4] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[3] is 0, but must be positive! (in '/tmp/RtmpAd2AVf/model-c2f548499100.stan', line 25, column 4 to column 38)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 1769.5 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 1860.2 seconds.
## Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 finished in 2844.4 seconds.
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 finished in 3053.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 2381.8 seconds.
## Total execution time: 3053.7 seconds.
## Warning: 6 of 2000 (0.0%) transitions hit the maximum treedepth limit of 10.
## See https://mc-stan.org/misc/warnings for details.
plot(precis( m_h2 , depth=2))
## 168 matrix parameters hidden. Use depth=3 to show them.
precis( m_h2 , depth=2)
## 168 matrix parameters hidden. Use depth=3 to show them.
## mean sd 5.5% 94.5% rhat ess_bulk
## bNl -0.7481725 0.08967661 -0.89094817 -0.6063615 1.002843 798.2844
## bNp 1.1878297 0.14180814 0.95764220 1.4083081 1.002752 1030.5424
## ap -2.3930813 0.13605202 -2.59388330 -2.1783546 1.007757 610.8038
## al 4.1783673 0.02957207 4.13421340 4.2240747 1.009969 565.9312
## sigma_id[1] 0.5163979 0.10217924 0.37901208 0.6919352 1.005957 713.7998
## sigma_id[2] 0.1205353 0.02323810 0.09008362 0.1606623 1.006809 701.4872
## sigma_id[3] 0.4916765 0.11387326 0.33270949 0.6959382 1.001755 1077.1723
## sigma_id[4] 0.3306998 0.07481651 0.23149843 0.4657552 1.001728 970.4360
precis( m_h2 , depth=3 , pars="Rho_id")
## mean sd 5.5% 94.5% rhat ess_bulk
## Rho_id[1,1] 1.000000000 0.0000000 1.0000000 1.00000000 NA NA
## Rho_id[2,1] -0.434492558 0.1901663 -0.7016454 -0.09527151 1.0020829 929.3918
## Rho_id[3,1] -0.420023675 0.2088108 -0.7197414 -0.06584689 0.9990519 1663.1131
## Rho_id[4,1] -0.191375372 0.2201639 -0.5342442 0.16532271 1.0038083 968.7778
## Rho_id[1,2] -0.434492558 0.1901663 -0.7016454 -0.09527151 1.0020829 929.3918
## Rho_id[2,2] 1.000000000 0.0000000 1.0000000 1.00000000 NA NA
## Rho_id[3,2] 0.252294131 0.2168762 -0.1204801 0.57163875 0.9998656 1875.1295
## Rho_id[4,2] 0.004964414 0.2124145 -0.3322853 0.33817939 1.0016527 1015.1019
## Rho_id[1,3] -0.420023675 0.2088108 -0.7197414 -0.06584689 0.9990519 1663.1131
## Rho_id[2,3] 0.252294131 0.2168762 -0.1204801 0.57163875 0.9998656 1875.1295
## Rho_id[3,3] 1.000000000 0.0000000 1.0000000 1.00000000 NA NA
## Rho_id[4,3] -0.029169033 0.2590645 -0.4565290 0.38869499 1.0026738 778.8528
## Rho_id[1,4] -0.191375372 0.2201639 -0.5342442 0.16532271 1.0038083 968.7778
## Rho_id[2,4] 0.004964414 0.2124145 -0.3322853 0.33817939 1.0016527 1015.1019
## Rho_id[3,4] -0.029169033 0.2590645 -0.4565290 0.38869499 1.0026738 778.8528
## Rho_id[4,4] 1.000000000 0.0000000 1.0000000 1.00000000 NA NA
Let plot the probability of being terrestrial in each context with got plots for each individual. Unsurprisingly, they are more terrestrial when using tools.
post <- extract.samples(m_h2)
dens(logistic(post$ap) , xlim=c(0.05, 0.5 ) , xlab="probability of being terrestrial")
dens(logistic(post$ap + post$bNp) , add=TRUE , col="slateblue" )
legend(x = "topright", legend = c("tool use", "non tool use") , fill = c("slateblue","black") , bty='n' )
blog=1.5
for(i in 1:17){
squanch <- logistic(post$ap + post$a_id[,i,1] )
points( median(squanch) , (i-.1)/blog , pch=18)
segments(x0=HPDI(squanch)[1] , x1=HPDI(squanch)[2] , y0=(i-.1)/blog , y1=(i-.1)/blog )
squanch <- logistic(post$ap + post$a_id[,i,1] + post$bNp + post$a_id[,i,2])
points( median(squanch) , (i+.1)/blog , col="slateblue" , pch=16)
segments(x0=HPDI(squanch)[1] , x1=HPDI(squanch)[2] , y0=(i+.1)/blog , y1=(i+.1)/blog , col="slateblue" )
}
Ok now we can plot the heights in decimeters (Sorry). Same story. Note 50 is 5 m.
dens( (1-logistic(post$ap))*exp(post$al) , xlim=c(0, 200) , xlab="height in dm during follows" , show.HPDI=.99999 )
dens((1-logistic(post$ap + post$bNp))*exp(post$al + post$bNl) , add=TRUE , col="slateblue" , show.HPDI=.99999)
legend(x = "topright", legend = c("tool use", "non tool use") , fill = c("slateblue","black") , bty='n' )
blog=200
for(i in 1:17){
squanch <- (1-logistic(post$ap + post$a_id[,i,1] ))*exp(post$al+ post$a_id[,i,3])
points( median(squanch) , (i-.1)/blog , pch=18)
segments(x0=HPDI(squanch)[1] , x1=HPDI(squanch)[2] , y0=(i-.1)/blog , y1=(i-.1)/blog )
squanch <- (1-logistic(post$ap + post$a_id[,i,1] + post$bNp + post$a_id[,i,2]))*exp(post$al + post$a_id[,i,3]+ post$bNl + post$a_id[,i,4])
points( median(squanch) , (i+.1)/blog , col="slateblue" , pch=16)
segments(x0=HPDI(squanch)[1] , x1=HPDI(squanch)[2] , y0=(i+.1)/blog , y1=(i+.1)/blog , col="slateblue" )
}
One of Tati’s big questions was, does the relationship between terrestriality out of tool use context, and terrestriality within tool use context differ? So if we are trying to predict terrestriality rates inside and outside of tool use context, are the the same within a quadrat.